Contents

Comparisons with PIOMAS sea ice thickness estimatesΒΆ

In this notebook, we compare our ICESat-2-derived sea ice thickness estimates with model-derived sea ice thickness from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS, Zhang and Rothrock, 2003)

import xarray as xr # For working with gridded climate data 
import pandas as pd
import numpy as np
import itertools
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting

# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection 
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore')

Read in dataΒΆ

Below, we’ll read in the data and restrict it to the Inner Arctic region. For the PIOMAS sea ice thickness, we’ll set 0 thickness to nan to match the ICESat2 data; the ICESat-2 dataset defines anywhere with no sea ice as nan, whereas PIOMAS sets anywhere with no sea ice to 0m thickness.

# Read/download the data 
book_ds = read_book_data() 

# Restrict to the Inner Arctic domain 
innerArctic = [1,2,3,4,5,6]
book_ds = book_ds.where(book_ds.region_mask.isin(innerArctic))

# Apply a filter to only seelct grid-cells where both ICESat-2 and PIOMAS have positive thickness (e.g. thickness greater than 1 cm)
book_ds = book_ds.where((book_ds["ice_thickness_int"]>0.01) & (book_ds["piomas_ice_thickness"]>0.01))

# Set 0 --> nan for PIOMAS thickness
#book_ds["piomas_ice_thickness"] = book_ds["piomas_ice_thickness"].where(book_ds["piomas_ice_thickness"]!=0, np.nan) 

# Apply an additional filter where the gridded monthly sea ice concentration is greater than 90% - i.e. just focus on the pack ice. 
book_ds = book_ds.where(book_ds["sea_ice_conc"]>0.9)

# Years over which to perform analysis
years = [2018,2019,2020] 

# Select variables of interest 
#book_ds = book_ds[["ice_thickness_int","ice_thickness_unc", "piomas_ice_thickness", "sea_ice_conc"]] # Grab data variables of interest

# Apply the thickness uncertainity as uper and lower values to simplify plotting uncertainity shading with Holoviews
book_ds['ice_thickness_upper']= book_ds.ice_thickness_int+book_ds.ice_thickness_unc
book_ds['ice_thickness_lower']= book_ds.ice_thickness_int-book_ds.ice_thickness_unc

Map monthly PIOMAS dataΒΆ

data_var = "piomas_ice_thickness"
#interactiveArcticMaps(book_ds[data_var], frame_width=500)

Map mean winter differencesΒΆ

Here, we’ll compute winter means for PIOMAS and ICESat-2 sea ice thickness. Then, we’ll use the interactiveArcticMaps function to generate maps of the winter means, along with a comparison map showing the gridcell differences.

is2_winter_means = compute_gridcell_winter_means(book_ds["ice_thickness_int"], years=years, start_month="Nov")
pio_winter_means = compute_gridcell_winter_means(book_ds["piomas_ice_thickness"], years=years, start_month="Nov")

for i in range(len(years)): 
    data1 = is2_winter_means.isel(time=i)
    data2 = pio_winter_means.isel(time=i)
    pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=data1.time.values.item())
    pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title=data2.time.values.item()) 
    diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="gridcell difference", title="Gridcell difference") 
    data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
    data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
    display(data_pl)

Mean monthly sea ice thickness comparisonΒΆ

Here, we’ll generate an interactive lineplot of the mean monthly sea ice thickness for PIOMAS and ICESat-2.

# Generate monthly means
mean_monthly_data = book_ds.mean(dim=("x","y"), keep_attrs=True)

# Holoviews requires us to generate individual line plots
is2_lineplot_p = mean_monthly_data.ice_thickness_int.hvplot.line(grid=True, label="ICESat-2", frame_width=500, frame_height=250) * mean_monthly_data.ice_thickness_int.hvplot.scatter(marker='s') # Overlay scatter plot to add markers
pio_lineplot_p = mean_monthly_data.piomas_ice_thickness.hvplot.line(grid=True, label="PIOMAS") * mean_monthly_data.piomas_ice_thickness.hvplot.scatter(marker='o') # Overlay scatter plot to add markers
# Generate shading based on the ICESat-2 uncertainty
is2_lineplot_p_unc = mean_monthly_data[["ice_thickness_lower", "ice_thickness_upper"]].hvplot.area(grid=True,y="ice_thickness_lower", y2 = "ice_thickness_upper", label="ICESat-2 uncertainty", alpha=0.1)

# Combine plots and display
comb_plot = (is2_lineplot_p*pio_lineplot_p*is2_lineplot_p_unc).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))
comb_plot.opts(title="PIOMAS and ICESat-2 sea ice thickness (Sep 2019 - Apr 2021)")
#comb_plot.opts(ylabel="effective mean thickness (meters)")
display(comb_plot)

A quick comparison between PIOMAS and ICESat-2 ice thickness shows that PIOMAS tends to have a lower mean thickness in the beginning of winter than ICESat2. Below, we’ll use the interactiveArcticMaps function to visualize any interesting patterns in the differences for September 2019, the month with the large difference between PIOMAS and ICESat-2 mean ice thickness (0.9 meters vs. 1.38 meters).

Note the large regional differences! Try playing around with different months and see how things look.

# Grab data for a given month
date = "Sep 2019"
pio = book_ds["piomas_ice_thickness"].sel(time=date)[0]
is2 = book_ds["ice_thickness_int"].sel(time=date)[0]

# Generate maps and display combined figure
pl_pio = interactiveArcticMaps(pio, vmin=0, vmax=4, colorbar=False, title="PIOMAS ice thickness ("+date+")") 
pl_is2 = interactiveArcticMaps(is2, vmin=0, vmax=4, title="ICESat-2 ice thickness ("+date+")") 
diff_pl = interactiveArcticMaps((pio - is2), vmin=-1.5, vmax=1.5, clabel="gridcell difference", cmap="coolwarm", title="Gridcell difference") 
data_pl = (pl_pio+pl_is2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
display(data_pl)